Path Integral Monte Carlo calculation of momentum distribution in solid 4 He 
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Abstract 

We perform calculations of the momentum distribution n(k) in solid 4 He by means of path integral Monte Carlo 
methods. We see that, in perfect crystals, n{k) does not depend on temperature T and that is different from the 
classical Gaussian shape of the Maxwell-Boltzmann distribution, even though these discrepancies decrease when the 
density of the system increases. In crystals presenting vacancies, we see that for T > 0.75 K, n(k) presents the same 
behavior as in the perfect crystal, but, at lower T, it presents a peak when k — > 0. 
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I. INTRODUCTION 



The observation, made by Kim and Chan in 2004, of a non classical rotational inertia (NCRI) in a 
torsional oscillator (TO) containing solid helium iQlQ] has generated, in the last years, a huge interest in the 
debate about supersolidity. The first theoretical speculations about the supersolid phase, i.e. a phase where 
crystalline order coexists with superfluidity and Bose-Einstein condensation (BEC), date back almost forty 
years y|45|]. Nevertheless, we are still far from a complete description of this phenomenon. The evidences 
of NCRI in solid helium have been confirmed in several TO experiments performed by other groups|@,0], 
which often gave controversial results. Nowadays, there is an overall agreement of all the data concerning 
the onset temperature of this phenomenon (T c ~ lOOmK), but the values of the superfluid density p. s /p 
reported so far can vary more than one order of magnitude, according to experimental conditions such as 
the way the crystal is prepared, its subsequent annealing or the 3 He concentration|8|]. These discrepancies 
suggest that the quality of the solid sample plays a very important role in these experiments. So, it is 
very important to understand the behavior of crystalline defects and to study if they are necessary to have 
superfluidity in a crystal and if they can be stable in the ground state of this system. 

From the theoretical point of view, the presence of BEC in a crystal can be detected by computing the 
momentum distribution w(k). The macroscopic occupation of the lowest energy state, in strongly interact- 
ing system like 4 He, appears in n(k) as a delta-peak for k = and a divergent behavior n(k) ~ \/k when 
k — > 0. Equivalently, one can get information about the BEC properties of a quantum system from the 
asymptotic behavior of the one-body density matrix pi(r), which is the inverse Fourier transform of n(k). 
If at large distances p\ (r) reaches a plateau the system present off-diagonal long range order (ODLRO), the 
condensate fraction no being this asymptotic value «o = lim r ^c»Pi(r). The theoretical study of solid 4 He 
at low temperatures cannot be developed analytically via a perturbative approach. It is therefore necessary 
the use of microscopic approaches to provide a reliable description of this phenomenon. Quantum Monte 
Carlo (QMC) methods have been extensively used for this purpose, but, so far, they have not been able 
to reproduce the experimental findings on supersolid 4 He. Path integral Monte Carlo simulations (PIMC) 
results have shown that a commensurate perfect crystal does not exhibit either superfluid fraction||9j] or 
ODLRO [ 10]. Instead, a non-zero condensate fraction has been observed in simulations of crystals present- 
ing defects, such as vacancies 111 ill or grain boundaries [12], and in simulations of an amorphous state of the 
solid jl3l] . Nevertheless, these simulations are not able to provide results in complete agreement with the 
experimental ones and therefore does not bring to a definitive answer of the supersolidity problem. 

In this paper, we perform PIMC calculations of the one-body density matrix and of the momentum 
distribution in solid 4 He by means of a very accurate sampling scheme. We will study commensurate and 
incommensurate hep crystals, both at zero and finite temperature. In section HI1 we describe the methods 
used in our simulations. In section [TTlJ we present the results obtained, at first, for perfect hep crystals at 
different densities, and subsequently, for hep crystals presenting vacancies. Finally, our conclusions are 
comprised in section HVl 



II. THE PATH INTEGRAL MONTE CARLO METHOD 

The path integral Monte Carlo (PIMC) method provides a fundamental approach in the study of strongly 



interacting quantum system at finite temperature 11411 . It is well known that the partition function 



Z = Tr(e-£ ft ) = / dR(R|e-' 3ft |R) (1) 



allows for a full microscopic description of the properties of a given system with Hamiltonian H = K + V 
at a temperature T = (^/J) -1 (the complete basis we use is the position basis \R) = |n 3 . . . r#) where the N 
particles are labeled). The noncommutativity of the kinetic and the potential energy operators (respectively, 
K and V) makes impractical a direct calculation of Z from Eq. Q] 
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The basic idea of PIMC is to use the convolution property of the thermal density matrix p(R,R';f}) = 
(R\e~P H \R'), in order to rewrite the partition function as 

M-l 

11 dRi p (R h R i+1 ; e) , (2) 
i=0 

with £ = j6/M and the boundary condition R M = Rq. For sufficiently large M, we recover the high- 
temperature limit for the thermal density matrix, where it is easy to separate the kinetic contribution from 
the potential one (Primitive Approximation). If one ignores the quantum statistics of the particles, the distri- 
bution law appearing in Eq. [2]is positive definite and can be interpreted as a probability distribution function 
which can be sampled by standard metropolis Monte Carlo methods. 

In practice, the PIMC method consists in mapping the finite-temperature quantum system to a classical 
system made up of closed ring polymers. This technique may be referred as an "exact" method, in the 
sense that using an accurate approximation for the high-temperature density matrices, the results are not 
affected by this approximation within the statistical error. However, its disadvantage is that the number M 
of convolution terms (beads) necessary to reach the convergence of Eq. |2]to the exact value of Z is inversely 
proportional to the temperature of the system: this means that, when approaching the interesting quantum 
regime at very low temperature, M increases fast making simulations hard, if not impossible, due to the very 
low efficiency in the sampling of the long chains involved. 

To overcome this problem, it is important to develop high-order approximation schemes for the density 
mauix, able to work with larger values of e. The approximation we use in this work is called Chin Approx- 



imation (CA)[15]. CA is based on a fourth order expansion of the e~$ H which makes use of the double 
commutator [[V',^],V], this term being related to the gradient of the interatomic potential. With respect to 
Takahashi-Imada Approximation!! 1 6ll . which is accurate to fourth order only for the trace, the new feature 
appearing in the CA is the presence of coefficients weighting the different terms in the expansion of the 
action: these coefficients are continuously tunable, making possible to force the error terms of fourth order 
to roughly cancel each other and get an effective sixth-order approximation. 

An additional problem we have to deal with when simulating quantum many-body systems with PIMC 
arises from the indistinguishable nature of the particles. If we deal with bosons like 4 He, the indistinguisha- 
bility of particles does not change the positivity of the probability distribution in Eq. [2] and the symmetry of 
p(R,R';fi) can be recovered via the direct sampling of permutations between the ring polymers represent- 
ing the quantum particles. To this purpose, a very efficient sampling is performed by the Worm Algorithm 
(WA)[17]: the basic idea of this technique is to work in an extended configuration space, given by the union 
of the ensemble Z, formed by the usual closed-ring configurations, and the ensemble G, which is made up 
of configurations where all the polymers but one are closed. Thanks to the presence of an open polymer, we 
are able to search the atoms involved in a permutation cycle by means of single particle updates, which do 
not suffer of a low acceptance rate and guarantee an efficient and ergodic sampling of the bosonic permuta- 
tions. We have to notice that the probability distribution used to sample the configurations in G is not equal 
to the one appearing in Eq. [2] and, therefore, these configurations cannot be used to calculate diagonal prop- 
erties, such as the energy or the superfluid density. However, the G-configurations can be used to compute 
off-diagonal observables such as the one-body density matrix pi(ri,r'j). Furthermore, the WA being able 
to sample both diagonal and off-diagonal configuration, is able to give an estimation of the normalization 
factor of pi(ri,ij). In this way, we are able to compute the properly normalized one-body density matrix 
and so to avoid the systematical uncertainties introduced by a posteriori normalization factor. 

The PIMC technique can be extended to zero temperature, in the so called path integral Ground State 
(PIGS) method 111811 . Indeed, we can see that the same imaginary-time evolution operator appearing in the 
definition of the thermal density matrix, can be used to project a trial wave function *Pt- to the exact ground- 
state wave function *Po according to the formula 

*P Q (R) = lim [dR , p(R,R , ;P) x ¥ T (R') ■ (3) 
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FIG. 1: The one-body density matrix p\ (r) (left) and the momentum distribution n(k) (right) for a commensurate hep 
crystal at density p = 0.0294 A -3 and at different temperatures: T = OK (circles), T = 1 K (squares) and T = 2K 
(diamonds). Statistical errors are below symbol size. 



As in the finite temperature approach, one can compute the ground-state averages of physical observables 
by factorizing p(R,R';(5) and studying the convergence with the number M of convolution terms. A good 
approximation for p(R,R';e) at small e, like CA, makes possible to choose = 1 in Eq. [3] and thus to 
obtain exact and completely model-independent results still with a very small number of beads M[19]. 



III. RESULTS 

We have carried out PIMC simulations of solid 4 He using a perfect hep lattice with a simulation box 
containing N = 180 atoms interacting with an Aziz pair potential[20]. At first, we perform calculations 
for a solid at a density p = 0.0294 A -3 , at three different temperatures. The results for Pi(r) and n{k) are 
shown in Figure [T] We can see that, at large r, p\ (r) decays exponentially, indicating that this system does 
not present ODLRO. We can also see that there is not a dependence of the momentum distribution with the 
temperature, in agreement with previous results by Clark and Ceperley floTl . 

In order to compare our results to the ones obtained by neutron scattering experiments, we compute 
the Compton profile of the longitudinal momentum distribution J(y), which is «(k) projected along the 
direction of the momentum Q of the incoming neutron and it is of easier experimental access than n(k). In 
the Impulse Approximation, which describe well the inelastic neutron scattering at high momentum transfer, 
J(y) and n(k) are related by the formula i2in 

J(y) = [ dkn(k)S(y-k Q ) =2% [ dkkn(k) (4) 

J ' J\y\ 

being k Q = k-^. 

In Figure|2l we compare the results for the longitudinal momentum distribution obtained from our n(k) at 
zero temperature with the fit obtained from experimental measurement by Diallo et al. for the same quantity 
in solid 4 He at molar volume V m = 20.01 cm 3 /mol (p = 0.0301 A~ 3 ) and a temperature T = 80nK[22]. We 
can see that our result are in a good agreement with the experimental ones. 

In Figure [2 we have also plotted the J(y) obtained by Ceperley in a PIMC simulation for a bec crystal 
at a density close to the melting of solid 4 He p = 0.0288A~ 3 and at a temperature T = 1.67K[23]. The 
difference between our momentum distribution and the one computed by Ceperley has to be attributed to the 
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FIG. 2: The longitudinal momentum distribution J(y): our results at T = OK and p — 0.0294 A 3 (circles) are com- 
pared with the fit obtained by Diallo et al ATM from neutron scattering experiment at T = 80nK and p = 0.0301 A -3 
(solid line) and with the PIMC results of CeperleyJH for a bcc solid at T = 1 .67 K and p = 0.0288 A~ 3 . 



larger density of the crystal we are simulating, and indicates that a larger coordination between the atoms 
in the solid cause a depletion of the low momentum states. 

This behavior is confirmed when simulating a crystal at higher densities. Figure [3] shows P\{r) and n(k) 
in a crystal at the density p = 0.0335 A -3 . For this density, the solid phase is stable over a larger range of 
temperature, making us able to simulate the system at temperatures up to T = 3 K. Comparing the results 
for the two different densities, we see that, when p increases, the one-body density matrix decays faster 
to zero and that the occupation of the low momentum states is appreciably decreased. It is also interesting 
to notice that p\(r) and n(k) are nearly independent of T even for temperatures larger than 7\, that is the 
superfluid transition temperature in the liquid phase. 

It is important to notice that at both densities, the shape of n(k) differs significantly from the classical 
Gaussian Maxwell-Boltzmann distribution. This feature is clearly shown in Figure [4] where we plot, on 
a logarithmic scale, n(k) as a function of k 2 . From this graph, it is easy to see the differences between 
the momentum distributions we obtained from PIMC and the straight line which represent a Gaussian 
distribution. We also notice that these differences are smaller in the system at p = 0.0335 A -3 , indicating 
that solid helium becomes more classic when the density increases. 

Finally, in order to give a deeper insight in the debate about supersolidity, it is important to understand 
how the momentum distribution changes when vacancies are present inside the crystal. For this purpose, 
we compute pi (r) and n(k) for a system made up of N = 179 atoms in a box which fits an N s = 180 sites 
hep lattice. We perform our calculations at the density p = 0.0294 A -3 and over a range of temperatures 
from r = 0Kte>r=2K. The results presented in figure [5] show that Pi(r) at T = OK presents a non-zero 
asymptote at large r, indicating that BEC is present inside the system. Our estimation of the condensate 
fraction is no = (9.0 ±0.8) x 10~ 4 , in agreement with Diffusion Monte Carlo results[24]. 

Instead, at finite temperature, we do not see a clear signal of ODLRO since all the pi (r) decay exponen- 
tially at large r: this indicates that the onset temperature T c at which the condensate appears in the system 
is below the lowest temperature we study, that is T c < 0.5 K. Nevertheless, from the study at finite temper- 
ature we notice an interesting dependence of P\{r) and n(k) with T: we can see that, for T > 0.75 K, Pi(r) 
differs only slightly from the same quantity calculated for the perfect crystal, while no difference can be 
seen for n(k) in the two systems. Instead, at T = 0.5 K, we see that the incommensurate crystal behaves dif- 
ferently since an additional peak centered in k = appears in the momentum distribution. This result may 
be explained assuming that at temperatures T > 0.75 K, the vacancy creates only a local distortion of the 
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FIG. 3: The one-body density matrix p\{r) (left) and the momentum distribution n(k) (right) for a commensurate 
hep crystal at density p = 0.0335 A -3 and at different temperatures: T = OK (circles), T = 1 K (squares), T = 2K 
(diamonds) and T = 3K (triangles). The dashed lines represent the same quantities computed for hep crystal at 
T = IK and p = 0.0294 A -3 




FIG. 4: The momentum distribution n(k) as a function of k 2 at T = 1 K and densities p = 0.0294 A 3 (solid line) and 
p = 0.0335 A -3 (dashed line). The dotted line represent a Gaussian n(k) and it is used to guide the eye. 



lattice which does not affect the momentum distribution, while, below a critical temperature Tq < 0.75 K, 
the defect begins to delocalize and to allow a larger occupation of the low momentum states. 



IV. CONCLUSIONS 

To summarize, we have computed the momentum distribution in solid 4 He by means of PIMC methods. 
These calculations are of fundamental importance in the study of BEC properties of quantum solids. The 
use of an effective sixth-order approximation for the action allows us to study the system by means of a 
very accurate sampling scheme, and the implementation of Worm Algorithm makes possible to calculate 
one-body density matrices which are correctly normalized. We have seen that our results are in good 
agreement with the ones obtained in neutron scattering experiments and indicates that solid 4 He is highly 



6 




FIG. 5: The one-body density matrix pi(V) (left) and the momentum distribution n(k) (right) for a hep crystal pre- 
senting a vacancy at density p = 0.0294A~ 3 and at different temperatures: T = OK (circles), T = 0.5 K (squares), 
T = 0.75 K (diamonds), T = IK (triangles up) and T = 2K (triangles down). The dashed lines represent the same 
quantities computed for a commensurate hep crystal at T = IK and p = 0.0294A -3 . Statistical errors are below 
symbol size. 



anharmonic, even though its behavior approaches the classical one at high densities. We have also shown 
that the presence of defects like vacancies affects the momentum distribution only at very low temperatures 
(T < 0.75 K). Calculations of pi(r) and n(k) in defected crystals at lower temperatures are now under way. 
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